Dimensiones del dataset, número de variables continuas y categóricas, distribuciones y comentarios generales
import pandas as pd
import numpy as np
pd.set_option('display.max_rows', None, 'display.max_columns', None)
fev = pd.read_csv(r'FEV_data.csv')
fev.head()
| seqnbr | subjid | age | fev | height | sex | smoke | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 301 | 9 | 1.708 | 57.0 | 2 | 2 |
| 1 | 2 | 451 | 8 | 1.724 | 67.5 | 2 | 2 |
| 2 | 3 | 501 | 7 | 1.720 | 54.5 | 2 | 2 |
| 3 | 4 | 642 | 9 | 1.558 | 53.0 | 1 | 2 |
| 4 | 5 | 901 | 9 | 1.895 | 57.0 | 1 | 2 |
fev.info()
fev.shape
<class 'pandas.core.frame.DataFrame'> RangeIndex: 654 entries, 0 to 653 Data columns (total 7 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 seqnbr 654 non-null int64 1 subjid 654 non-null int64 2 age 654 non-null int64 3 fev 654 non-null float64 4 height 654 non-null float64 5 sex 654 non-null int64 6 smoke 654 non-null int64 dtypes: float64(2), int64(5) memory usage: 35.9 KB
(654, 7)
Después de cargar el dataset y ver la información del mismo, vemos que tenemos 654 observaciones y 7 columnas, sin embargo la primera (seqnbr) es unicamente el número de secuencia o index de cada observación, mientras que subjid es el ID del paciente (en la buena teoría debería ser un número random). Por otro lado las columnas de sex y smoke son categóricas (realmente binarias), por lo que vamos a convertirlas al tipo 'category', no sin antes asegurarnos de que efectivamente contienen unicamente dos valores diferentes cada una.
fev.nunique()
seqnbr 654 subjid 654 age 17 fev 575 height 56 sex 2 smoke 2 dtype: int64
factors = list(fev.loc[:, fev.nunique()<= 10])
fev[factors] = fev[factors].astype('category')
fev.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 654 entries, 0 to 653 Data columns (total 7 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 seqnbr 654 non-null int64 1 subjid 654 non-null int64 2 age 654 non-null int64 3 fev 654 non-null float64 4 height 654 non-null float64 5 sex 654 non-null category 6 smoke 654 non-null category dtypes: category(2), float64(2), int64(3) memory usage: 27.2 KB
Con esto, tendríamos una variable objetivo (fev) y 4 son posibles variables predictoras, de las cuales tenemos 2 variables continuas (age y heigth) y 2 variables categóricas.
No tenemos missing values puesto que todas las columnas tienen 654 datos, sin embargo debemos comprobar que no haya outliers, además queremos visualizar la distribución de las variables. Para esto haremos una tabla con datos descriptivos y posteriormente una inspección gráfica.
fev.describe()
| seqnbr | subjid | age | fev | height | |
|---|---|---|---|---|---|
| count | 654.00000 | 654.000000 | 654.000000 | 654.000000 | 654.000000 |
| mean | 327.50000 | 37169.571865 | 9.931193 | 2.636780 | 61.143578 |
| std | 188.93782 | 23690.860350 | 2.953935 | 0.867059 | 5.703513 |
| min | 1.00000 | 201.000000 | 3.000000 | 0.791000 | 46.000000 |
| 25% | 164.25000 | 15811.000000 | 8.000000 | 1.981000 | 57.000000 |
| 50% | 327.50000 | 36071.000000 | 10.000000 | 2.547500 | 61.500000 |
| 75% | 490.75000 | 53638.500000 | 12.000000 | 3.118500 | 65.500000 |
| max | 654.00000 | 90001.000000 | 19.000000 | 5.793000 | 74.000000 |
fev.describe(exclude = np.number)
| sex | smoke | |
|---|---|---|
| count | 654 | 654 |
| unique | 2 | 2 |
| top | 1 | 2 |
| freq | 336 | 589 |
import seaborn as sns
import matplotlib.pyplot as plt
import plotly.express as px
def histogram_boxplot(data, xlabel = None, title = None, font_scale=2, figsize=(9,8), bins = None):
# Crear ventana para los subgráficos
f2, (ax_box2, ax_hist2) = plt.subplots(2, sharex=True, gridspec_kw={"height_ratios": (.15, .85)}, figsize=figsize)
# Crear boxplot
sns.boxplot(x=data, ax=ax_box2)
# Crear histograma
sns.histplot(x=data, ax=ax_hist2, bins=bins) if bins else sns.histplot(x=data, ax=ax_hist2)
# Pintar una línea con la media
ax_hist2.axvline(np.mean(data),color='g',linestyle='-')
# Pintar una línea con la mediana
ax_hist2.axvline(np.median(data),color='y',linestyle='--')
# Asignar título y nombre de eje si tal
if xlabel: ax_hist2.set(xlabel=xlabel)
if title: ax_box2.set(title=title, xlabel="")
# Mostrar gráfico
plt.show()
def plot(col):
if col.dtypes != 'category':
print('Variable Continua')
histogram_boxplot(col, xlabel = col.name, title = 'Distibución continua')
else:
print('Cat')
fig = px.bar(col.value_counts())
fig.show()
fev.drop(['seqnbr', 'subjid'], axis = 1).apply(plot)
Variable Continua
Variable Continua
Variable Continua
Cat
Cat
age None fev None height None sex None smoke None dtype: object
No tenemos missing data, no vemos errores en los datos; parece haber algunos outliers en la variable fev, sin embargo al tratarse de la variable objetivo y además, al no estar tampoco tan alejados de los demás valores al punto de parecer irreal, seguiremos trabajando con ellos para no perder posible información valiosa.
En cuanto a la distribución, tenemos para age y height distribuciones normales, fev, no es tan simétrica, pero de igual forma, posiblemente esté muy cerca de la normalidad. Mientas que las variables categóricas, tenemos por un lado Sex, que está muy bien distribuida entre sus dos categorías, con cantidad de datos muy similar para cada una; y por el otro lado, smoke (que como sería de esperarse al tratarse de niños y adolescentes) tiene muy poca representación de pacientes fumandores, la mayoría está en la categoría 2 (asumiremos que es no fumadores).
Parece ser una base de datos bastante limpia y lista para modelar!
Como se mencionó antes, la variable seqnbr, es el número de secuencia, por lo que no trabajaremos con ella.
La variable subjid pareciera ser también algun número de identificación del paciente, en ese caso no debería tomarse en cuenta puesto que se trata de números aleatorios.
La variable smoke, está muy poco representada en una de sus dos categorías, por lo que podría aportarnos muy poca información e incluso información sesgada, no vamos a descartarla desde ya, sino que trabajaremos con ella para ver si de alguna forma nos aporta información valiosa, pero está en la mira para descartarla en caso de que no nos aporte mucho.
Vamos a hacer una inspección visual y general, para hacernos una idea de las posibles relaciones entre las variables continuas
g = sns.PairGrid(fev.iloc[:, 2:])
g.map_upper(sns.scatterplot)
g.map_diag(sns.histplot)
g.map_lower(sns.regplot)
<seaborn.axisgrid.PairGrid at 0x1c818948250>
Parece haber una clara relación entre las variables de edad y altura con la variable objetivo y también parece haber una dependencia entre ellas; sin embargo aún no hemos visto nada con respecto a sex y smoke que son categóricas.
Utilicemos pandas profiling para tener más información:
import pandas_profiling
pandas_profiling.ProfileReport(fev.iloc[:, 2:])
Summarize dataset: 0%| | 0/5 [00:00<?, ?it/s]
Generate report structure: 0%| | 0/1 [00:00<?, ?it/s]
Render HTML: 0%| | 0/1 [00:00<?, ?it/s]
En este reporte vemos resumido mucho de lo que analizamos antes, sin embargo hay datos adicionales en la sección de correlaciones: por un lado nos confirma una correlación fuerte entre age y heigth con la variable objetivo (y entre ellas también!) y por otro lado, nos indica que hay cierta relación entre sex y smoke con fev, aunque no es tan fuerte. Para esta segundo punto (correlación categórica-numérica) el reporte utiliza el método v de cramer; por lo que sería interesante ver este análisis en detalle.
Primero vamos a separar la variable objetivo para poder manipular las demás variables sin riesgo de modificar o transformar la variable de respuesta
varObj = fev.fev
inputs = fev.drop(['seqnbr', 'fev', 'subjid'],axis=1)
inputs.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 654 entries, 0 to 653 Data columns (total 4 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 age 654 non-null int64 1 height 654 non-null float64 2 sex 654 non-null category 3 smoke 654 non-null category dtypes: category(2), float64(1), int64(1) memory usage: 11.9 KB
Y ahora si, usaremos la función de VCramer para un análisis preliminar de la relación de las variables con la variable objetivo (aunque ya tenemos una ligera idea de esto con base en la inspección visual que hicimos, pero queremos hacer la comparación de la correlación para cada input)
import scipy.stats as stats
def cramers_v(var1, vObj):
if not var1.dtypes == 'category':
bins = min(5,var1.value_counts().count())
var1 = pd.cut(var1, bins = 5)
if not vObj.dtypes == 'category':
bins = min(5,vObj.value_counts().count())
vObj = pd.cut(vObj, bins = 5)
data = pd.crosstab(var1, vObj).values
vCramer = stats.contingency.association(data, method = 'cramer')
return vCramer
# Ejemplo uso univariante
#cramers_v(vinosCompra['Etiqueta'],vinosCompra['Beneficio'])
# Aplicar la función al input completo contra la objetivo
tablaCramer = pd.DataFrame(inputs.apply(lambda x: cramers_v(x,varObj)),columns=['VCramer'])
# Obtener el gráfico de importancia de las variables frente a la objetivo continua según vcramer
import plotly.express as px
px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones respecto al FEV').update_yaxes(categoryorder="total ascending")
Ahora si, podemos ver ordenados los inputs según su correlación con la variable objetivo!
A lo que vinimos, vamos a generar el modelo de regresión lineal para la predicción del FEV; lo primero es separar los datos en dos grupos: training y test; posteriormente generamos el modelo con los datos de training
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(inputs, varObj, test_size=0.2, random_state=42)
# Comprobamos dimensiones
print('Training dataset shape:', X_train.shape, y_train.shape)
print('Testing dataset shape:', X_test.shape, y_test.shape)
Training dataset shape: (523, 4) (523,) Testing dataset shape: (131, 4) (131,)
# Unimos los datos de training de la variable objetivo, con los datos de training de los inputs:
fev_training = X_train.join(y_train)
fev_training.head()
| age | height | sex | smoke | fev | |
|---|---|---|---|---|---|
| 321 | 12 | 71.0 | 1 | 2 | 4.550 |
| 456 | 12 | 68.0 | 1 | 2 | 4.411 |
| 340 | 10 | 62.0 | 1 | 2 | 1.937 |
| 29 | 9 | 60.0 | 2 | 2 | 2.100 |
| 570 | 10 | 65.0 | 1 | 2 | 3.090 |
# Generamos el modelo completo:
from statsmodels.formula.api import ols
reg_lineal = ols('fev ~ age + height + sex + smoke', data=fev_training).fit()
reg_lineal.summary()
| Dep. Variable: | fev | R-squared: | 0.766 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.764 |
| Method: | Least Squares | F-statistic: | 423.4 |
| Date: | Sun, 05 Mar 2023 | Prob (F-statistic): | 1.09e-161 |
| Time: | 09:45:19 | Log-Likelihood: | -281.13 |
| No. Observations: | 523 | AIC: | 572.3 |
| Df Residuals: | 518 | BIC: | 593.6 |
| Df Model: | 4 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -4.3576 | 0.271 | -16.067 | 0.000 | -4.890 | -3.825 |
| sex[T.2] | -0.1511 | 0.038 | -4.028 | 0.000 | -0.225 | -0.077 |
| smoke[T.2] | 0.0652 | 0.066 | 0.983 | 0.326 | -0.065 | 0.196 |
| age | 0.0632 | 0.011 | 5.678 | 0.000 | 0.041 | 0.085 |
| height | 0.1041 | 0.005 | 19.150 | 0.000 | 0.093 | 0.115 |
| Omnibus: | 11.617 | Durbin-Watson: | 2.092 |
|---|---|---|---|
| Prob(Omnibus): | 0.003 | Jarque-Bera (JB): | 17.038 |
| Skew: | 0.178 | Prob(JB): | 0.000200 |
| Kurtosis: | 3.810 | Cond. No. | 932. |
Como habíamos mencionado antes, la variable smoke está muy poco representada en una de sus categorías (es normal, los adolescentes no deberían fumar), además vemos que el p-value de su coeficiente de correlación es mayor a 0.05 (es el único mayor a 0); por lo que, como habíamos indicado antes, vamos a excluirla del modelo y ver cual es el efecto en el R cuadrado:
reg_lineal1 = ols('fev ~ age + height + sex', data=fev_training).fit()
reg_lineal1.summary()
| Dep. Variable: | fev | R-squared: | 0.765 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.764 |
| Method: | Least Squares | F-statistic: | 564.2 |
| Date: | Sun, 05 Mar 2023 | Prob (F-statistic): | 6.86e-163 |
| Time: | 09:45:52 | Log-Likelihood: | -281.61 |
| No. Observations: | 523 | AIC: | 571.2 |
| Df Residuals: | 519 | BIC: | 588.3 |
| Df Model: | 3 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -4.2829 | 0.260 | -16.450 | 0.000 | -4.794 | -3.771 |
| sex[T.2] | -0.1540 | 0.037 | -4.117 | 0.000 | -0.227 | -0.081 |
| age | 0.0601 | 0.011 | 5.629 | 0.000 | 0.039 | 0.081 |
| height | 0.1043 | 0.005 | 19.223 | 0.000 | 0.094 | 0.115 |
| Omnibus: | 11.573 | Durbin-Watson: | 2.100 |
|---|---|---|---|
| Prob(Omnibus): | 0.003 | Jarque-Bera (JB): | 17.840 |
| Skew: | 0.155 | Prob(JB): | 0.000134 |
| Kurtosis: | 3.850 | Cond. No. | 892. |
El R cuadrado del modelo practicamente no cambió al excluir smoke, por lo que lo mantendremos así. Parece ser que tenemos un buen modelo con un R2 de 0.765; sin embargo solo por entenderlo un poco mejor, vamos a ver cual es el aporte de cada variable al R2:
import patsy
# Generamos las matrices de diseño según la fórmula de modelo completo
y, X = patsy.dmatrices('fev ~ age + height + sex + smoke', fev_training, return_type='dataframe')
from relativeImp import relativeImp
names=X.columns.tolist()[1:]
# Calculamos importancia relativa
df_results = relativeImp(X.join(y), outcomeName = 'fev', driverNames = names)
# Ordenamos valores
df_results.sort_values(by='normRelaImpt', ascending=False)
| driver | rawRelaImpt | normRelaImpt | |
|---|---|---|---|
| 3 | height | 0.443669 | 57.937958 |
| 2 | age | 0.277787 | 36.275651 |
| 0 | sex[T.2] | 0.024235 | 3.164859 |
| 1 | smoke[T.2] | 0.020075 | 2.621533 |
px.bar(df_results,x='normRelaImpt',y='driver',title='Importancia relativa por aportación al R2').update_yaxes(categoryorder="total ascending").show()
Según esto, podríamos eliminar (además de smoke) también a sex del modelo, veamos como cambia:
reg_lineal2 = ols('fev ~ age + height', data=fev_training).fit()
reg_lineal2.summary()
| Dep. Variable: | fev | R-squared: | 0.758 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.757 |
| Method: | Least Squares | F-statistic: | 812.9 |
| Date: | Sun, 05 Mar 2023 | Prob (F-statistic): | 8.87e-161 |
| Time: | 09:46:25 | Log-Likelihood: | -290.02 |
| No. Observations: | 523 | AIC: | 586.0 |
| Df Residuals: | 520 | BIC: | 598.8 |
| Df Model: | 2 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -4.5960 | 0.253 | -18.181 | 0.000 | -5.093 | -4.099 |
| age | 0.0523 | 0.011 | 4.900 | 0.000 | 0.031 | 0.073 |
| height | 0.1095 | 0.005 | 20.417 | 0.000 | 0.099 | 0.120 |
| Omnibus: | 18.099 | Durbin-Watson: | 2.110 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 31.130 |
| Skew: | 0.231 | Prob(JB): | 1.74e-07 |
| Kurtosis: | 4.103 | Cond. No. | 852. |
El R2 se disminuye ligeramente, pero nada considerable, por lo que tenemos dos posibles modelos lineales: uno que incluye edad, altura y sexo, y otro únicamente con edad y altura (que desde hace rato vimos que eran las que parecían tener mayor correlación).
Vamos a evaluar ambos modelos con los datos de test:
from sklearn.metrics import mean_squared_error, r2_score
# Predicciones para test modelo completo
fev_pred1 = reg_lineal1.predict(X_test)
# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred1)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred1))
Mean squared error: 0.40 Coefficient of determination: 0.80
# Predicciones para test modelo completo
fev_pred2 = reg_lineal2.predict(X_test)
# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred2)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred2))
Mean squared error: 0.41 Coefficient of determination: 0.79
Tenemos datos muy similares para ambos modelos; los dos son aceptables, pero una variable menos nos simplificaría mucho el modelo. Veamos las interacciones:
Ya habíamos visto en nuestro análisis gráfico, que hay una aparente interacción entre age y height, que además como vimos antes, vienen siendo las variables más relevantes para la predicción del valor de FEV, por lo que vamos a tomar en cuenta esa interacción:
reg_lin_interac2 = ols('fev ~ age + height + age*height', data= fev_training).fit()
reg_lin_interac2.summary()
| Dep. Variable: | fev | R-squared: | 0.783 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.781 |
| Method: | Least Squares | F-statistic: | 622.7 |
| Date: | Sun, 05 Mar 2023 | Prob (F-statistic): | 1.72e-171 |
| Time: | 09:48:36 | Log-Likelihood: | -261.65 |
| No. Observations: | 523 | AIC: | 531.3 |
| Df Residuals: | 519 | BIC: | 548.3 |
| Df Model: | 3 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -0.4138 | 0.593 | -0.698 | 0.485 | -1.579 | 0.751 |
| age | -0.4306 | 0.063 | -6.789 | 0.000 | -0.555 | -0.306 |
| height | 0.0410 | 0.010 | 4.013 | 0.000 | 0.021 | 0.061 |
| age:height | 0.0077 | 0.001 | 7.712 | 0.000 | 0.006 | 0.010 |
| Omnibus: | 17.104 | Durbin-Watson: | 2.071 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 37.858 |
| Skew: | -0.053 | Prob(JB): | 6.02e-09 |
| Kurtosis: | 4.314 | Cond. No. | 2.25e+04 |
reg_lin_interac1 = ols('fev ~ age + height + sex + age*height', data= fev_training).fit()
reg_lin_interac1.summary()
| Dep. Variable: | fev | R-squared: | 0.785 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.784 |
| Method: | Least Squares | F-statistic: | 474.1 |
| Date: | Sun, 05 Mar 2023 | Prob (F-statistic): | 1.49e-171 |
| Time: | 09:48:40 | Log-Likelihood: | -258.17 |
| No. Observations: | 523 | AIC: | 526.3 |
| Df Residuals: | 518 | BIC: | 547.6 |
| Df Model: | 4 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -0.5429 | 0.592 | -0.918 | 0.359 | -1.705 | 0.619 |
| sex[T.2] | -0.0968 | 0.037 | -2.635 | 0.009 | -0.169 | -0.025 |
| age | -0.3881 | 0.065 | -5.961 | 0.000 | -0.516 | -0.260 |
| height | 0.0431 | 0.010 | 4.229 | 0.000 | 0.023 | 0.063 |
| age:height | 0.0071 | 0.001 | 6.971 | 0.000 | 0.005 | 0.009 |
| Omnibus: | 14.277 | Durbin-Watson: | 2.066 |
|---|---|---|---|
| Prob(Omnibus): | 0.001 | Jarque-Bera (JB): | 28.474 |
| Skew: | -0.054 | Prob(JB): | 6.56e-07 |
| Kurtosis: | 4.138 | Cond. No. | 2.26e+04 |
De nuevo vemos que la presencia o no de sex en el modelo, no afecta mucho, de hecho afecta aún menos en el modelo que incluye la interacción de edad y altura (que por cierto, parece ser un modelo levemente mejor que el que no la incluye si nos basamos en el R2.
Vamos a ver la evaluación de estos modelos al comparar las predicciones con los datos de test y ya luego veremos la validación cruzada de ambos modelos:
# Predicciones para test modelo completo
fev_pred3 = reg_lin_interac1.predict(X_test)
# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred3)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred3))
Mean squared error: 0.39 Coefficient of determination: 0.82
# Predicciones para test modelo completo
fev_pred4 = reg_lin_interac2.predict(X_test)
# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred4)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred4))
Mean squared error: 0.39 Coefficient of determination: 0.81
Efectivamente vemos que el aporte de sex al modelo, es mínimo, podemos excluirlo para simplificar el modelo sin perder información importante!
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from sklearn.linear_model import LinearRegression
def cross_val_lin(formula, data, seed=12345):
y, X = patsy.dmatrices(formula, data, return_type='dataframe')
model = LinearRegression()
# Establecemos esquema de validación fijando random_state (reproducibilidad)
cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=seed)
# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X, y, cv=cv)
# Sesgo y varianza
print('Modelo: ' + formula)
print('Coeficiente de determinación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
#sns.violinplot(y=scores,palette='viridis')
return(scores)
form1 = 'fev ~ age + height'
form2 = 'fev ~ age + height + age*height'
# Creamos lista de fórmulas
list_form = [form1,form2]
list_form
# Aplicamos a toda la lista la función creada (devuelve un dataframe pero está transpuesto)
list_res = pd.DataFrame(map(lambda x: cross_val_lin(x,fev, seed=2022),list_form))
list_res
Modelo: fev ~ age + height Coeficiente de determinación R2: 0.760 (0.030) Modelo: fev ~ age + height + age*height Coeficiente de determinación R2: 0.783 (0.033)
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | 91 | 92 | 93 | 94 | 95 | 96 | 97 | 98 | 99 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.787667 | 0.716968 | 0.799950 | 0.729255 | 0.761944 | 0.791089 | 0.776350 | 0.688730 | 0.788862 | 0.765495 | 0.772354 | 0.786487 | 0.769036 | 0.718342 | 0.762030 | 0.784979 | 0.741495 | 0.785038 | 0.709679 | 0.772243 | 0.750218 | 0.789058 | 0.750301 | 0.784801 | 0.719166 | 0.711838 | 0.761917 | 0.725076 | 0.818966 | 0.767828 | 0.723102 | 0.712695 | 0.774206 | 0.809178 | 0.785401 | 0.823126 | 0.691026 | 0.811060 | 0.708003 | 0.747786 | 0.797421 | 0.736606 | 0.727761 | 0.768323 | 0.780752 | 0.777619 | 0.757366 | 0.727513 | 0.764900 | 0.772503 | 0.802393 | 0.762426 | 0.778155 | 0.698572 | 0.755787 | 0.803356 | 0.710792 | 0.782174 | 0.763558 | 0.727841 | 0.809943 | 0.768817 | 0.756976 | 0.748240 | 0.746915 | 0.73496 | 0.757857 | 0.791671 | 0.743574 | 0.768408 | 0.794476 | 0.765307 | 0.778449 | 0.736414 | 0.724770 | 0.760133 | 0.738737 | 0.783343 | 0.806688 | 0.720035 | 0.726006 | 0.721742 | 0.775253 | 0.781166 | 0.775853 | 0.731595 | 0.778746 | 0.72327 | 0.784665 | 0.798909 | 0.771248 | 0.750733 | 0.771886 | 0.727637 | 0.769702 | 0.775609 | 0.761466 | 0.763527 | 0.750511 | 0.760841 |
| 1 | 0.808889 | 0.743159 | 0.815954 | 0.749891 | 0.798461 | 0.811031 | 0.812616 | 0.721601 | 0.811489 | 0.762009 | 0.801723 | 0.817343 | 0.795215 | 0.740051 | 0.777243 | 0.813767 | 0.770766 | 0.816728 | 0.697378 | 0.800250 | 0.794001 | 0.801966 | 0.779349 | 0.807909 | 0.730690 | 0.716065 | 0.801887 | 0.733507 | 0.845703 | 0.796791 | 0.737427 | 0.746988 | 0.790016 | 0.819267 | 0.825149 | 0.840542 | 0.706053 | 0.840481 | 0.726434 | 0.780371 | 0.834036 | 0.756675 | 0.741305 | 0.791445 | 0.800239 | 0.813321 | 0.780701 | 0.724106 | 0.796349 | 0.794709 | 0.810594 | 0.769931 | 0.805100 | 0.734141 | 0.790558 | 0.836799 | 0.744945 | 0.803868 | 0.771287 | 0.752825 | 0.822430 | 0.779777 | 0.786041 | 0.788226 | 0.765887 | 0.75224 | 0.770672 | 0.815152 | 0.788136 | 0.789405 | 0.821415 | 0.789580 | 0.812760 | 0.758850 | 0.733953 | 0.779442 | 0.765022 | 0.805122 | 0.824413 | 0.753467 | 0.731746 | 0.748743 | 0.796179 | 0.817650 | 0.799745 | 0.743433 | 0.785447 | 0.77693 | 0.823026 | 0.790566 | 0.796992 | 0.751288 | 0.806779 | 0.759025 | 0.791933 | 0.797006 | 0.793718 | 0.756039 | 0.793153 | 0.785192 |
results = list_res.T.melt()
results.columns = ['Modelo','R2']
results.groupby(['Modelo']).describe().T
| Modelo | 0 | 1 | |
|---|---|---|---|
| R2 | count | 100.000000 | 100.000000 |
| mean | 0.760046 | 0.782857 | |
| std | 0.030248 | 0.032920 | |
| min | 0.688730 | 0.697378 | |
| 25% | 0.736050 | 0.756516 | |
| 50% | 0.764229 | 0.790287 | |
| 75% | 0.781418 | 0.807061 | |
| max | 0.823126 | 0.845703 |
Vemos un mejor coeficiente de determinación para el modelo con interacciones, aunque también una leve desmejora en el error estandar, pero esta desmejora es mínima, entonces nos vamos a dejar influenciar más por el coeficiente de determinación (que tampoco tiene diferencia abismal, pero si se ve mejoría). Veamos de una forma gráfica la diferencia en el R2:
sns.boxplot(x='Modelo',y='R2',data=results,palette='viridis')
<AxesSubplot:xlabel='Modelo', ylabel='R2'>
Por lo que podemos ver, ambos modelos son relativamaente similares, aunque el modelo con form2 (con interacción) es el de mejor ajuste por lo que este sería nuestro modelo final!
form2 = 'fev ~ age + height + age*height'